工具&方法 | 6行代码教你用python做OLS回归（内附CFPS实例） | 自由微信

查看原文

其他

工具&方法 | 6行代码教你用python做OLS回归（内附CFPS实例）

Original 企研数据 ·王学习数据Seminar 2021-06-04

收录于话题

#启研学社 157

#软件应用 104

写在前面的话

阅读本文前需要掌握的基础知识：Python 的基础知识、 numpy 的基础知识、 pandas 的基础知识、基本的计量知识。如果你还不会，那么本文也会介绍一些 python 语法的基础内容，方便大家理解。

随着数据资源的日渐丰富，学者们越来越多的需要接触到大数据的处理，许多学者还是习惯使用 Stata 对数据进行处理，而 Stata 由于其自身的限制，在处理大数据集时要么需要强劲的处理性能（昂贵的硬件成本），要么需要等待较长时间（更加昂贵的时间成本）。Python 和 R 也就日渐进入学者的视野，相对于 R ，Python 的语法更为简单，成为一部分学者的首选。

在数据处理上，numpy 和 pandas 的组合，使得 Python 能够轻松应对千万级别的数据处理。在攻克数据处理这一环后，在数据应用上，除了新潮的机器学习、深度学习的方法，对于现阶段社科学者来说，计量可能才是最现实的。在 Python 中处理的数据如果还需要调回到 Stata 中做计量，那未免太「蹩脚」。今天，数据Seminar 公众号将带大家体验 Python 上的第三方计量库：Statsmodels[1] 。

简介

Statsmodels 是一个Python的第三方模块，他封装了许多计量模型，方便学者直接调用。所谓封装，就相当于 Stata 中一个 reg 命令，代表了最基础的 OLS 回归命令，在Statsmodels 中也有类似 reg 的语句，提供给 OLS 估计。另外 Statsmodels 的开源协议为 BSD [2]（基本上对于用户来说属于为所欲为协议，你可以任意使用这款扩张包，具体参见链接地址）。

statsmodels is a Python module that provides classes and functions for the estimation of many different statistical models, as well as for conducting statistical tests, and statistical data exploration. An extensive list of result statistics are available for each estimator. The results are tested against existing statistical packages to ensure that they are correct. The package is released under the open source Modified BSD (3-clause) license. The online documentation is hosted at statsmodels.org.

安装

如果你使用 Anaconda 安装的 python ，那么：

conda install statsmodels

如果你使用 pip 管理你的python包，那么：

pip install statsmodels # python2

或者：

pip3 install statsmodels #python3

“精读”代码，初探：OLS

官方的说明文档[3]：http://www.statsmodels.org/stable/gettingstarted.html

首先先看了一下 OLS 的示例：

import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf

dat = sm.datasets.get_rdataset("Guerry", "HistData").data
results = smf.ols('Lottery ~ Literacy + np.log(Pop1831)', data=dat).fit()
print(results.summary())

运行后：

本文采用的编辑器为 Google 的 colab

是的，OLS就是这么几行，具体来看内容，其实进行 OLS 回归只需要一行内容。为了方便 Python 基础不那么好的同学，我们驻行“精读”一遍这些代码：

1. 导入模块

import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf

首先是需要导入 numpy、statsmodels 的相关模块，numpy 推荐的缩写命名 np，statsmodels 推荐的缩写命名是 sm，这样后面的代码就可以通过 np.[具体的函数\方法] ，sm. [具体的函数\方法] 的方式来更加简洁的调用了。

2. 生成测试数据集

dat = sm.datasets.get_rdataset("Guerry", "HistData").data

由于并没有现成数据，我们可以通过 statsmodels 中自带的数据集，关于该数据集的具体信息可以参考[4]https://vincentarelbundock.github.io/Rdatasets/doc/HistData/Guerry.html，简单来说，这是一个关于犯罪、识字等相关内容的社会科学领域的数据集。

这个数据集以DataFrame 格式的储存在了 dat 这个变量中。DataFrame 格式在 Python 中被广泛的使用，可以说无论是机器学习还是本文介绍的计量相关内容，DataFrame 格式相关的操作都需要熟练掌握。关于 DataFrame 的相关内容，需要更加细致的去学习Python的另一个第三方库 pandas ，本文不再赘述。

3. 回归操作

results = smf.ols('Lottery ~ Literacy + np.log(Pop1831)', data=dat).fit()

首先是smf.ols， smf 就是前面的 statsmodels.formula.api ，OLS 回归的方法在 statsmodels.formula.api下，如果前文没有声明 smf 代表statsmodels.formula.api 的话，你这里可就要写成 statsmodels.formula.api.ols 了，可见合理的缩写声明对于代码的简洁美观非常有帮助。

接下来是括号内的内容('Lottery ~ Literacy + np.log(Pop1831)', data=dat)这其实是OLS模型需要传入的参数，以,为区隔，第一个参数是回归模型，第二个参数是数据。

在这行代码中，回归模型是：'Lottery ~ Literacy + np.log(Pop1831)'，数据是前面加载的dat ，比较方便的是，statsmodels，直接从数据dat中读取了变量名称，并不需要进行额外的赋值操作。我们具体看一下这个模型的写法，对于没有接触过 Python 和 R 的读者来说，这种写法可能相对有些陌生，其实这个模型代表了：

这种写法参考了R的相关内容，具体是由 Python 中的 Patsy 库实现的。上面这个就是一个计量模型（他当然不是计算 Literacy + np.log(Pop1831)的和，Literacy 和 np.log(Pop1831) 都是矩阵）。这就和 Stata 中 reg y x0 x1 x2 是一样，只是在Python 中模型需要写成y~x0+x1+x2，而 reg 需要写成 smf.ols()。

这个模型中还用到了一个知识点：np.log ，就是说引用了 numpy 中的 log 函数，对变量Pop1831 取对数，就是使用 numpy 进行 log 运算。如果你直接写 log(Pop1831)，Python是不知道 log 是做什么的，所以要告诉 Python，这个 log 来自于 numpy，这样 Python 就能正常处理 Pop1831 这个数据切片（“切片” 是Python数据处理中常用名词，可以理解为一个 DataFrame 的一列、一行或者几列、几行的数据）了。

最后 data=dat 就是声明对于这个模型，使用的数据时 dat，这样前面的模型中的 Lottery 、Literacy 和 Pop1831，statsmodels 都能智能的去从 DataFrame 切片（切片操作需要学习Pandas哦），然后应用到模型中。

括号内的内容就是进行回归操作的核心。

括号外还有.fit()，就是告诉 Python，可以进行回归计算了。如果没有这个.fit() ，Python 只会记录这个模型和相关的数据信息，而添加了.fit()，记录时回归后的结果。

上图可以看到储存在计算机中的内容是有差异的

同样，.fit() 后输出的结果也并非像 Stata 那样那样直接输出的是结果，而是将结果储存，可以进行多种方式的调用。

如果想输入结果，那么就使用 results.summary() ，将其打印出来，正如示例中展示的那样：

OLS 实例

这里简单演示一下 OLS 回归，数据处理这一部分内容，在 Python 中主要使用 numpy 和 pandas，这里就不演示，本文直接读取处理完的数据来演示 OLS 的内容：

本文的数据来自CFPS，主要研究的内容是子女认知能力和父母是否创业的关系，处理后的第一个回归模型的内容主要是，母亲的创业状态（self_employed_m_16），教育年限(edu)，是否参加课外补习(extra_classes)，家庭药品支出(expend_medical)、自评健康状况(health)、家庭氛围(quarrel)、与母亲同住时长(livewith_m)、与父亲同住时常(livewith_f)、log的家庭收入(lg_fin)、子女数量(num_chd)、母亲的受教育程度(edu_m)、父亲的受教育程度(edu_f)、户口(hukou)、东中西部虚拟变量（east 、 west），然后被解释变量就是认知能力(y_rznl)。

import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf

df = pd.read_csv("testing.csv")
OLS_1 = smf.ols("y_rznl~ self_employed_m_16 + edu + extra_classes + \
expend_medical + health + quarrel + livewith_m + livewith_f + \
lg_fin + num_chd + edu_m + edu_f + hukou + east + west", data=df).fit()
print(OLS_1.summary())

最后输出结果：

而stata输出的结果（置信区间是97.5%）：

可以看到，结果基本一致。

最后，在OLS这里案例中，可以看到 Python 输出结果基本与 Stata 相同，在代码的撰写上， Python 同样非常方便，并且由于 Python 是一门编程语言，在一些循环的撰写上非常有优势，例如我需要多次循环某几个变量，以及对应的被解释变量：

import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf

df = pd.read_csv("testing.csv")
y_list =[y1, y2, y3, y4] # 被解释变量的列表
x0_list = [x0_1, x0_2, x0_3, x0_4, x0_5] # 某一个解释变量的列表
model_list = []
for y in y_list:
for x0 in x0_list:
model = "%s ~ %s + edu + extra_classes + \
expend_medical + health + quarrel + livewith_m + livewith_f + \
lg_fin + num_chd + edu_m + edu_f + hukou + east + west" % (y, x0)
OLS = smf.ols(model, data = df).fit()
print(OLS.summary())